Bold diagrammatic Monte Carlo: 
A generic technique for polaron (and many-body?) problems 
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We develop a Monte Carlo scheme for sampling series of Feynman diagrams for the proper self- 
energy which are self-consistently expressed in terms of renormalized particle propagators. This 
approach is used to solve the problem of a single spin-down fermion resonantly interacting with 
the Fermi gas of spin-up particles. Though the original series based on bare propagators are sign- 
alternating and divergent one can still determine the answer behind them by using two strategies 
(separately or together): (i) using proper series re-summation techniques, and (ii) introducing renor- 
malized propagators which are defined in terms of the simulated proper self-energy, i.e. making the 
entire scheme self-consistent. Our solution is important for understanding the phase diagram and 
properties of the BCS-BEC crossover in the strongly imbalanced regime. On the technical side, we 
develop a generic sign-problem tolerant method for exact numerical solution of polaron- type models, 
and, possibly, of the interacting many-body Hamiltonians. 

PACS numbers: 05.30.Fk, 05.10.Ln, 02.70.Ss 



I. INTRODUCTION 

Modern science has radically changed our view of the 
physical vacuum. Instead of naive "empty space" we 
have to deal with a complex groundstate of an interact- 
ing system, and, strictly speaking, there is no fundamen- 
tal difference between the outer space, helium, or any 
other condensed matter system. With this point of view 
comes understanding that the notion of a "bare" particle 
is somewhat abstract since its coupling to vacuum fluc- 
tuations, or environment, may strongly (sometimes radi- 
cally) change particle properties at energies addressed by 
the experimental probes. The polaron problem [l| is by 
now canonical across all of physics with the same ques- 
tions about particle energy, effective mass, residue, etc., 
being asked for different types of particles, environments 
and coupling between them In a broader context, 
particles are not necessarily external objects unrelated 
to a given vacuum, but quasiparticle excitations of the 
very same ground state. Thus, the solution of the po- 
laron problem paves the way to the effective low energy 
theory of a given system, and, to large extent, determines 
basic properties of all condensed matter systems at low 
temperature. 

At the moment, there is no generic analytic or numeric 
technique to study quasiparticle properties for arbitrary 
strongly interacting system. Analytic solutions are typi- 
cally (with few exceptions, see, e.g., Q) based on pertur- 
bative corrections to certain limiting cases 0, H, [j, [H, [f| 
or variational treatment 0- Several numeric schemes 
were suggested in the past, but all of them have limita- 
tions either in the system size, system dimension, inter- 
action or environment type. Exact diagonalization and 
variational methods in the low-energy subspace 0, Q are 
mostly restricted to one-dimensional models with short- 
range interactions. The continuous time path integral 
formulation [lfj| works for the lattice model with linear 



coupling between the particle and bosonic environment, 
but it can not be generalized to fermionic environment, 
sign-alternating coupling (as in the t-J model [Hi]), nor 
is it suitable for continuous-space models. 

In this article (which follows a short communication 
[T3|). we develop a Monte Carlo technique which sim- 
ulates series of Feynman diagrams for the proper self- 
energy. The diagrammatic Monte Carlo (Diag-MC) tech- 
nique was used previously to study electron-phonon po- 
larons [HI, [ijj] • The essence of Diag-MC is in interpret- 
ing the sum of all Feynman diagrams as an ensemble 
averaging procedure over the corresponding configura- 
tion space. It was considered essential for the method to 
work that the series of diagrams for the Green's function 
be convergent and sign-positive. Though the configura- 
tion space of diagrams for polarons in Fermi systems is 
more complex, similar methods of generating and sam- 
pling the corresponding configuration space can be used. 
The crucial difference is that in the Fermi system we have 
to deal with the sign-alternating and divergent (at least 
for strong coupling) series. Under these conditions a di- 
rect summation of all relevant Feynman diagrams for the 
Green's functions is not possible, and one has to develop 
additional tools for (i) reducing the number of diagrams 
by calculating self-energies rather than Green's functions, 
(ii) employing the "bold-line" trick in the form of the 
Dyson equation which allows self-consistent summation 
of infinite geometric series and further reduces the num- 
ber of self-energy diagrams, and, if necessary, (iii) extrap- 
olating Diag-MC results to the infinite diagram order for 
a divergent series. At the moment we do not see any ob- 
vious limitations of the new method since even divergent 
sign alternating series can be dealt with to obtain reliable 
results. We believe that our findings are important in a 
much broader context since the Diag-MC approach to the 
many-body problem has essentially the same structure. 

As a practical application of the method we consider a 
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particle coupled to the ideal Fermi sea via a short-range 
potential with divergent s-wave scattering length. This 
problem was recognized as the key one for understanding 
the phase diag ram of the population imbalanced resonant 
Fermi gas [l5l [l6|. In particular, to construct the energy 
functional describing dilute solutions of minority (spin 
down) fermions resonantly coupled to the majority (spin 
up) fermions one has to know precisely the quasiparticle 
parameters of spin-down fermions since they determine 
the coefficients in the energy expansion in the spin-down 
concentration x^: the linear term is controlled by the 
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polaron energy, and the x/ term is determined by the 
polaron mass [16j |. 

The general Hamiltonian we deal with in this article 
can be written as 

H = Hp - A B /2m + J drV{r - R) n(r) , (1) 

where Hp is the Hamiltonian of the ideal spin-up Fermi 
gas with density n and Fermi momentum kp, R is the 
particle coordinate, and V(r — R) is the interaction po- 
tential of finite range r between the particle and the 
spin-up Fermi gas. In what follows we refer to (TT]) as the 
Fermi-polaron problem. The specifics of the BCS-BEC 
crossover physics in the strongly imbalanced regime is 
two-fold: one is that the particle and the Fermi gas have 
the same bare mass m (in what follows we will use units 
such that m = 1/2 and kp = 1), and the other is that 
one has to take explicitly the so-called zero-range reso- 
nant limit when ro — > 0, but the s-scattering length a 
remains finite, i.e. kpa remains fixed for kp ro — > 0. In 
this limit, the nature of the interaction potential is irrel- 
evant, and the same results will be obtained e.g. for the 
neutron matter and Cesium atoms. We note, however, 
that the method we develop for the numeric solution of 
the resonant Fermi-polaron problem is absolutely gen- 
eral and can be used for an arbitrary model described by 
Eq. Q). 

It turns out that the structure of the phase diagram 
is very sensitive to polaron parameters. If the state with 
a dilute gas of spin-down fermions is stable at all val- 
ues of kpa then the solution of the single-particle prob- 
lem would define the phase diagram in the vicinity of 
the multicritical point discussed recently by Sachdev and 
Yang [It]], where four different phases meet. At this point 
the spin-down fermion forms a bound state with a spin- 
up fermion thus becoming a spin-zero composite boson 
("molecule"); i.e. quasi-particles radically change their 
statistics. The multicritical point, however, may be ther- 
modynamically unstable if the effective scattering length 
between the composite bosons and spin-up electrons is 
large enough, and the analysis of Refs. [HI, based on 
the fixed-node Monte Carlo simulations finds evidence in 
favor of this scenario. Phase separation was also found 
in calculations based on mean-field and narrow-resonance 
approaches (both at finite and zero temperature) see e.g. 
Refs. [H El, [H, SI El El, though with strong quan- 
titative deviations from results based on the fixed-node 



Monte Carlo simulations [19(. On the experimental side, 
MIT experiments [26[ are in good agreement with the 
predictions made in Ref. [l9[ , while Rice experiments [27| 
are not. The origin of discrepancy between the two ex- 
periments is not understood. Our results for polaron 
energies are in excellent agreement with Ref. [l9j]. 

It is worth noting that the model (TT]) (in general, the 
particle mass is different from that of the Fermi gas) is 
also known as the Anderson orthogonality problem with 
recoil [H, S3- It can be also considered as a specific 
example of a particle coupled to the Ohmic dissipative 
bath (see monograph [30| for numerous other examples 
and connections to realistic systems). 

The paper is organized as follows. In Sec. [TT] we dis- 
cuss the configuration space of Feynman diagrams for 
self-energy in momcntum-imaginary-time representation 
(both in the particle and molecule channels), and ex- 
plain how polaron parameters can be obtained in this 
representation. In Sec. IIIII we describe a Monte Carlo 
algorithm for generating and sampling the correspond- 
ing diagrammatic space. A small technical Section IIVI 
deals with numerically evaluating the effective T-matrix 
by bold diagrammatic Monte Carlo. We present and dis- 
cuss results in Sec.[V] In particular, we show that one can 
use re-summation techniques for divergent series of dia- 
grams based on bare propagators to determine the final 
answer. In Sec. I VII we further advance the algorithm by 
employing bold-line approach in which the entire scheme 
is self-consistently based on renormalized ("bold-line") 
propagators. We present our conclusions and perspec- 
tives for future work in Sec. IVIII 



II. CONFIGURATION SPACE OF 
SELF-ENERGY DIAGRAMS 

As mentioned above, when coupling between spin- 
down and spin-up fermions is strong enough they from 
a composite boson, or molecule state. In what follows, 
we will be using the term "polaron" in a narrow sense, 
i.e. only for the unbound fermionic spin-down excita- 
tion. For the composite boson we will be using the 
term "molecule". Since our goal is to calculate parti- 
cle properties for arbitrary coupling strength we have to 
consider one- and two-particle channels on equal footing. 
In the rest of this section we render standard diagram- 
matic rules for irreducible self-energy diagrams in both 
channels, with an emphasis on specifics of working in the 
imaginary-time representation. 



A. Polaron channel 

We start from the definition of the single-particle 
Green's function (see e.g. [3ll ]) 

G(T,r)=-<T T #r,r$(0)), (2) 
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and its frequency-momentum representation 



G(£,p)= / e l (« r - p r ) G(t, r) dvdr . (3) 

Here ip(r,r) is the fermion annihilation operator at the 
space-time point (r, r) . For the ideal spin- up Fermi gas 
at T = we have 



G T (£,p) = - 



1 



if — p 2 /2m + ep 



(4) 



The vacuum Green's function for the spin-down po- 
laron is 



g| 0) (t,P,m) 



-9{t) e 



-(p /2m— fi)j 



(•5) 



where 9 is the step function, and /i is a free parameter. 
From Dyson's equation for the polaron, see Fig. [TJ one 
finds 



Gi(£,p,m) 



^-p 2 /2m + M-S(f,p,/i) 



(6) 



where self-energy S is given by the sum of all irreducible 
diagrams (i.e. diagrams which can not be made discon- 
nected but cutting through the line) taken with 
the negative sign. Taking into account that in the r- 
representation both G[ and E depend on \x only through 
exponential factors exp(/zr), one obtains G|(f,p,/i) = 
- if*, P) and E(f , p, /x) = E(f - ifi, p). 
If polaron is a well-defined quasi-particle, then its en- 
ergy E(p) and residue Z(p) can be extracted from the 
asymptotic decay 



G 4 (r,p^)-f-Ze-( B -^, 



(7) 



This asymptotic behavior immediately implies that the 
function Gj_(f — ijJt, p) has a pole singularity 



- i/J,p) 



Z(p) 



regular part . (8) 



ii + H-E{p) 

and comparing the result 



Now setting /i = E(p) in 
to ©, we conclude that 



if /Z = if - p 2 /2m + J5 - E(0, p, E) + i£A(p, E) , (9) 
where (we take into account that dTi/d^ = idYl/dfj,) 

5E(0,p,/i) 



A(p,E) 



(10) 



fi=E 



This yields two important relations (see also [3l|): 

E = p 2 /2m + E(0, p, £) , (11) 

and 

1 



1 + A(p,£0 



(12) 



O 



-G 



FIG. 1: Dyson equation for the single-particle Green's func- 
tion. 



Equation (fTTj) allows one to solve for E provided 
E(t, p, /i) is known. All we have to do is to calculate 
the the zero-frequency value of E for fi = E 



/•oo 

E = p 2 /2m+ / E(r,p,^)e(' B -^ T rfT 



(13) 



After E is found, the residue is obtained from Eq. |T 
using 



/>oo 

A(p,E)=- rE(r,p, M )e( £ -^dr. 
Jo 



(14) 



Note also that the dependence on /i drops out from both 

(USD and mm. 

A comment is in order here. Strictly speaking, po- 
laron and molecule poles exist only for p = because 
the fermionic bath they couple to is gapless. However, 
the spectrum E(jp) is well-defined in the limit p — > 0, as 
the decay width vanishes faster than [E(p) — E(Q)]. To 
have stable quasi-particlcs, one can use a trick of intro- 
ducing a gap A in the environment spectrum, e.g. by 
redefining the dispersion relation for spin- up fermions: 
£k — > max(£k,A). In the p — > limit, the system- 
atic error vanishes faster than [E(p) — E(0)], provided 
A ~ [-E(p) — E(0)]. It should be also possible to work 
with A's essentially larger than [E(p) — E(0)] and extrap- 
olate to A — > 0. In particular, such an extrapolation is 
possible (and is implicitly implied) at the analytical level 
in the relation for the effective mass, which we consider 
below. 

One way to determine the effective mass is to calculate 
the quasi-particle energy as a function of momentum for 
a set of small but finite values of p and fit it with the 
parabola. It is, however, possible to skip this procedure 
and to write a direct estimator for the effective mass in 
terms of the calculated self-energy. Acting with the op- 
erator Vp on both sides of Eq. (fTTj) and taking the limit 
p — * we get 



l + A 1 



Bn 



(15) 



Bo = 3 J o dr e^-^ [V P E(r, p, M )]| p=Q , (16) 
where A = A(p = 0) and E = E(p = 0). 



B. Molecule channel 



In this case we start with the two-particle propagator 
A-(r,p) = -(T T $ p (T)*t(Q)), (17) 
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FIG. 2: Defining functions Q and K in the two-particle chan- 
nel. 



where 



rZq 
2^ 



VqV'T(p-q)V'i(q) , 



(18) 



and ifq is the pair wavefunction (in momentum represen- 
tation) that localizes the relative distance between two 
particles. If there is a bound state (molecule), then 



K(t,p,(i) -> -Z mo ie 



-(E mo i-n)T 



(19) 



and the pair propagator in the frequency representation 
has a pole: 



(p) 



i£ + n - E mo i(p) 



regular part . (20) 



Now we introduce yet another function, that features 
the same molecular pole, but has a simpler diagrammatic 
structure. The specifics of the resonant zero-range limit 
is that the sum of all ladder diagrams for the interac- 
tion potential V{r) has to be considered as a separate 
diagrammatic element. We denote this sum by T(r,p) 
and consider it as a "bare pair propagator" . Of course, 
the same approach can be taken in a general case to re- 
place the bare interaction potential with the scattering T- 
matrix, but in the zero-range limit we really do not have 
any other alternative for the ultra-violet-divergence-free 
formulation. The sum of ladder diagrams takes the ultra- 
violet physics into account exactly and allows to express 
r(r, p) in terms of the s-scattering length a. The ladder 
structure of diagrams absorbed in T(t, p) also explains 
why we treat it as a "pair propagator" (we will depict 
it with a double-line, see Fig [5]). The exact expression 
for r is readily obtained from the geometrical series (in 
frequency domain) 



-U 



-U) 2 U 



where the polarization operator is defined by 

dq/(2^) 3 



U(rj,p) 



q >k F q 2 /2m + (p - q) 2 /2m - 77 



(21) 



(22) 



and rj = uj + sf + (i. Using suitable ultra-violet regular- 
ization, one can cast the same expression in the universal 
form which depends only on the s-scattering length: 



r-\ v ,p) 



Aita 87r 



\J p 2 — 4rn7/ — 11(77, p) i (23) 



dq/(27r) 



q <k F q 2 /2m + (p - q) 2 /2m - 77 



(24) 



For finite density of spin-up fermions converting Eq. 
to the imaginary time domain has to be done numeri- 
cally. One possibility is to use the inverse Laplace trans- 
form. We employed the bold diagrammatic Monte Carlo 
technology [32] to achieve this goal and further details 
are provided in Sec. IIVI The two-dimensional function 
r(r, p) is tabulated prior to the polaron simulation. 

In Fig. [5] we define function Q that can be viewed as a 
renormalized pair propagator related to T by the Dyson 
equation. In the upper panel, we show the diagrammatic 
structure for K, which includes dotted lines represent- 
ing external functions <pq, grey squares representing sums 
of all T-irreducible diagrams, and the renormalized pair 
propagator Q. By T-irreducible diagrams we understand 
diagrams which can not be made disconnected by cutting 
them through a single T-line. All T-reducible diagrams 
are absorbed in the Q function which is shown in the 
lower panel in Fig. [2j The grey circle has nearly the 
same structure as the grey square (the zeroth order term 
is present in the crossed square, but not in the crossed 
circle): since T is defined as the sum of ladder diagrams, 
all terms which include ladder-type structures based on 
free one-particle propagators have to be excluded from Q 



and K. With the replacements — + Q, G 



(0) 



r, and 



£ — * K we find an exact analogy between the one- and 
two-particle propagators. 

The analogy can be carried out further by noting that 
the structure of diagrams in Fig. [2] implies that Q has the 
same pole as K, while the rest of the functions simply 
change the value of the quasi-particle residue. Thus, if 
molecule is a well-defined excitation we expect that 



Q(£-i/x,p) 



Z mo \(p) 



i£ + (J,- E mo i(p) 



regular part . (25) 



This explains why introducing the function Q is conve- 
nient: now Eqs. (fT2|) -(fT6 |) are immediately generalized to 
the molecule case (up to replacements mentioned above). 



III. WORM ALGORITHM FOR FEYNMAN 
DIAGRAMS 



In this Section we describe how the configuration space 
of Feynman diagrams for E and K is parameterized and 
updated using Diag-MC rules. Our algorithm is designed 
to treat polaron and molecule channels on equal footing. 
We achieve this goal by introducing auxiliary diagrams 
which contain two "loose" ends called "worms" — this was 
proven to be an efficient strategy for reducing the MC 
autocorrelation time when simulations are performed in 
the configuration space with complex topology [H, ■ 
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FIG. 3: The backbone of the cyclic diagram. 




FIG. 4: Forward connection. The arc represents — Gf(r = 

T c +T d + T e ). 



A. Cyclic diagrams 

We start with the introduction of cyclic diagrams. 
Though we work in the imaginary time representa- 
tion at T = when r 6 [0, oo), it is still conve- 
nient not to specify the time origin and to consider di- 
agrams on the imaginary time circle. The backbone 
of each cyclic diagram is a pre-diagram illustrated in 
Fig. [3l It consists of a cyclic chain of the struc- 
ture Gf ) (r a )r(T b )G; o) (r c )r(r d )G; o) (T e )r(r / ) ... (all 
the times are positive). We do not explicitly show "direc- 
tions" of the propagators, since these are unambiguously 
fixed by the global direction of all the backbone lines, 
which we select — without loss of generality — to be from 
right to left. With this convention, the left (free) spin- 
up end of any T-line is outgoing, while the right end is 
incoming. A physical diagram is obtained by pairwise re- 
placing free spin-up ends with propagators G j . There are 
two ways to connect incoming and outgoing lines: (i) for- 
ward, i.e., in the direction of the backbone propagators, 
and (ii) backward (opposite to forward). Forward (back- 
ward) connections result in propagators Gj with positive 
(negative) times, see Figs. 0] and [5] for illustrations. They 
represent particle (hole) excitations in the fermionic en- 
vironment. It is important to emphasize that in cyclic 
diagrams the only time-variables are the positive time- 
lengths of G { ° ] 's and T's. There is no absolute time, and, 
correspondingly, all moments in time are equivalent. 



B. Worms 

To have a MC scheme which simulates diagrams in 
one- and two-particle channels on equal footing we ex- 
tend the space of physical diagrams by allowing diagrams 
with two special end-points, or "worms" . They will be 
denoted X and M. and stand for the unconnected incom- 
ing (outgoing) spin-up ends, see Fig.[5]for an illustration. 
Correspondingly, the entire updating scheme is based on 



FIG. 5: Backward connection. The pair of F-ends being con- 
nected is precisely the same as in Fig. 3] but the direction is 
opposite, and the arc represents — G-\{j = — t/ — r a — Th)- 




FIG. 6: A diagram with two worms, I and M. 

manipulations with X and M.. As it will become clear 
soon, of special importance for normalization of MC re- 
sults is the first-order diagram with the worm, see Fig. [7] 
Its weight consists of just two factors, G|°^(r a ) and r(r&). 



C. Parametrization of diagrams 

Apart from the diagram order and its topology, we se- 
lect time intervals of T's and G^ 's, and momenta of the 
spin-up propagators as independent variables. The mo- 
menta of T's and 's are then unambiguously defined 
by the momentum conservation law, while the time in- 
terval of a spin-up propagator is obtained by summing 
up the time intervals of T's and G^'s covered by this 
propagator. Technically, we find it convenient to work ex- 
plicitly in the particle-hole representation for the spin-up 
propagators when backward spin- up propagator is under- 
stood as a forward hole propagator with the opposite mo- 
mentum. This is achieved by introducing a non-negative 
function, 

G(r,p) = (- G / T ' p )\ (26) 
I G T (-r,-p) , P < Pf , v ; 

which is assigned to all spin-up lines (the global fermionic 
sign of the diagram is defined separately, by standard 
diagrammatic rules). All momenta assigned to the spin- 
up lines are understood as momenta of the corresponding 
G-propagators; i.e. they are either momenta of particles 
(for forward propagators they are non-zero only for p > 
Pf), or momenta of holes (for backward propagators they 
are non-zero only for p < Pf). An explicit formula for G 
(subject to the above conditions) is 

G(r,p) = 6»(r)e- |p2/2m - Ei?|T . (27) 
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as a parameter. The acceptance ratio for this update is 



FIG. 7: "Normalization" diagram. It is the simplest diagram 
with the worm; its weight is a product of G| '(r a ) and r(rfe). 



To simplify the description of updates below we will 
generically refer to G , l°' ) - and T-propagators as backbone 
lines (BBLs) and denote them as D. The diagram order 
N is given by the total number of spin-down propagators. 
We also use special counters to characterize the topology 
of the diagram. For each BBL we define a cover num- 
ber ; N c , equal to the total number of G-lines covering a 
given BBL. A backbone line with N c = is called un- 
covered. Finally, physical diagrams — the ones without 
worms — are divided into polaron (0) and molecule (1) 
sectors; the diagram sector is defined by the difference 
between the number of particle and hole spin-up propa- 
gators covering any of the -lines (the same result is 
obtained by analyzing propagators covering T-lines after 
adding unity for the spin-up particle participating in the 
ladder diagrams). 



D. Updates 

The cyclic structure of diagrams in combination with 
the possibility of considering non-physical diagrams al- 
lows one to construct a very simple ergodic set of updates. 
The minimal set consists of two complementary pairs, In- 
sert/Delete and Open/Close, and one self-complementary 
update Reconnect. The description that follows aims at 
the most transparent and straightforward realization of 
updates. Standard performance enhancement tricks and 
optimization protocols are not mentioned. In particular, 
we base our considerations on the updating scheme which 
may propose a change leading to a forbidden configura- 
tion. Such proposals are rejected as if they result in zero 
acceptance ratio. 

Insert: This update applies only to physical — no 
worms — diagrams and is rejected otherwise. First, con- 
sider the polaron sector. Select at random one of the 
propagators; if it is covered, the update is automatically 
rejected. If the selected propagator is uncovered, insert 
a pair of new propagators, r(n,p) and G^(t 2 , p), right 
after the selected one. The new r(n, p)-propagator is 
supposed to contain X and M. at its ends. Worms radi- 
cally simplify this diagram-order increasing update since 
due to conservation laws the momenta of new BBLs are 
equal to the global momentum of the diagram p. The 
times n and t 2 are drawn from normalized probabil- 
ity distributions Wt{t\) and W|(t 2 ) (arbitrary at this 
point). Note that Wt{t\) and Wj(T 2 ) can depend on p 



(28) 



where Cat is an artificial weighing factor ascribed to all 
worm diagrams of order N (it can be used for optimiza- 
tion purposes and depend on p too). A natural choice 
for VF-functions is to make them proportional to BBL, 
i.e. 



C<"V.P) 



p)dT> 



Then, to have an acceptance ratio of order unity and 
independent of p we choose 

C N = j^, A = J r(r 1>P )dri J Gf ] (r 2 ,p)dr 2 . (30) 

In the rest of the paper, we will refer to this choice of 
W and Gat as "optimized", though we do not mean that 
it is the best one possible for the entire scheme. For the 
"optimized" choice 



P ins = N/(N+l) . 



(31) 



In the molecule sector, we essentially repeat all steps, 
up to minor modifications. Now the propagator being 
selected is T (once again, the update is rejected if the se- 
lected propagator is covered.) Then, a pair of new prop- 
agators, G^(ri,p) and r(r 2 ,p), is inserted in front of 
the selected uncovered propagator. The new propagator 
r(r 2 , p) inherits the outgoing spin-up line previously con- 
nected to the selected T; the latter gets instead the M.- 
end of the worm while the X-end is attached to the new 
r. The acceptance ratio is identical to (j28|) [for the opti- 
mized choice it is N/(N+ 1)]. The polaron and molecule 
sectors are mutually exclusive due to particle conserva- 
tion, thus only one type of the Insert update is applicable 
to a given diagram. 

Delete: This update converts worm diagrams to phys- 
ical ones while reducing the diagram order. It applies 
only to diagrams of order N > 1 with worms being sep- 
arated by one uncovered BBL. Otherwise, the update is 
rejected. If worms are separated by an uncovered BBL, 
then the left and right neighbor BBLs are also uncovered, 
and an update opposite to Insert is possible. In Delete 
we simply remove two consecutive BBL and worms from 
the diagram. Its acceptance probability is the inverse of 

Eq. m, 



Pad = 



1 



{N-1)C 



W r (Ti)W;(73) 



r(r 1 , P )G; o) (r 2 , P ) 



(32) 



With Eqs. Q29J), Q3UJ), we have 
.Pdd = N/(N - 1) . 



(optimized) (33) 
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Close: The update applies only to diagrams with the 
worms. The proposal is to connect X and M with a line 
G(t, q) and eliminate worms from the diagram. The mo- 
mentum variable q is drawn from the probability distri- 
bution W\ (q) , while r is the time interval between X and 
M. to be covered by the new propagator. There are two 
ways of connecting X and M, forwards and backwards. 
The ambiguity is automatically resolved by the absolute 
value of the momentum variable q: if q > pp (q < Pf), 
the propagator is supposed to go forwards (backwards). 
In practice, we first select the direction (with equal prob- 
abilities), and then generate the momentum variable q 
accordingly: cither q > pp or q < pp. 

The acceptance ratio for this update is 



G(r,q) jj A^r^p^) 



(2TrfNC N W^q) LI D v (tv,p„) ' 



(34) 



where the subscript v runs through all BBLs to be cov- 
ered by the new propagator (clearly, r = ~Y^ V t v ). Primes 
indicate new values of the corresponding momenta: 



(35) 



As usual, the distribution function Wj (q) can depend on 
t and the direction of the propagator. The natural choice 
for this function would be 



W T (q) 



G(r,q) 
Q(r) 



n(r) = 



I q > PF ( T > q) dc l (forward) , 
I q<PF G 1 ( T > 0) dc i (backward) 



leading to the optimized acceptance ratio 

.Pel 



2Afi(T)j-r£> 1/ (r„ > p / J 



(2tt) 3 



D u (t v ,p u ) 



(36) 



(37) 



(38) 



In this Section we deal with diagrams based on bare 
propagators. To avoid double-counting, we have to ex- 
clude all cases which can be reduced to ladders already 
summed in T's. When M. and X are on the nearest- 
neighbor BBL the proposal to connect them with the 
spin-up particle propagator is rejected. The last rule to 
be monitored it to restrict all physical diagrams to be 
either in the polaron or molecule sectors, i.e. sectors dif- 
ferent from and 1 are not allowed. 

Open: The update applies only to physical diagrams 
and proposes to create a worm by selecting at random 
and removing one of the spin-up propagators. The ac- 
ceptance ratio is given by the inverse of ([34]) . ([38]) 



= (2ir) 3 NC N Wf(q) jt DAt v , P ' v ) 
° P 2 G T (r,q) \ L D„(T v ,p v ) ' 



(39) 



where the subscript v runs through all BBLs covered by 
the propagator, r = ~Y^ v t v , ano - 1 ^ s the momentum of 



the selected spin-up propagator, 
values of the BBL momenta: 

pi = Vv + q 

In the optimized version, we have 



Primes indicate new 



(40) 



P — 

1 op — 



(2tt) 3 jt D v (t v ,p' v ) 
2AO(t) 11 D v { Tv ,p„) 



(41) 



Reconnect: The purpose of this update is to change 
the topology of diagrams with the worm. The proposal is 
to select at random one of the G-propagators and swap its 
outgoing end with M. ; the momentum of the propagator 
remains the same, only its time variable changes from To 
to Tq. The acceptance ratio is given by: 



G(r ,q) 11 A^.P") 



(42) 



The subscript v runs through all BBLs that will change 
their momenta (and cover numbers -/V c 's) as a result of 
the update. Topologically, there are two different situa- 
tions (the two a complementary to each other in terms 
of the update): (i) M. is covered by the propagator in 
question, and (ii) M. is not covered by the propagator. 
Correspondingly, the new values of the diagram variables 
are calculated as 



(t - t, p v + q) (i) , 
(t +T,p u - q) (ii) , 



(43) 



with t = J2v T v 

The above set of updates is ergodic. It can be sup- 
plemented by additional updates that may improve the 
algorithm performance by more efficient sampling of the 
diagram variables and topologies. Over-complete sets of 
updates are also useful for meaningful tests of the de- 
tailed balance. The possibilities are endless, and here we 
simply mention two updates we have been using. 

Time shift: We propose new time variables, t v — > r£, 
for all uncovered BBLs (labeled here with the subscript 
v). The acceptance probability is 



n 



W v (t v ,p) P) 
W v (ji,p) D v {t v ,p) 



(44) 



All uncovered propagators have the same momentum p. 
With the optimized choice for W u {t v , p) oc D u (t„, p) the 
acceptance ratio becomes unity. 

Redirect: Here we propose to select at random one 
of the G-propagators and change its direction to the op- 
posite. Simultaneously, we change the momentum of the 
selected propagator (resulting in new momenta for all 
BBL it covers or will cover as a result of the update) . Let 
the selected propagator be G(t, q), and the new one be 
G(t', q') with t = t„, and t' = ^ A t\, where v runs 
through all BBLs covered by the propagator Gt (t, q) and 
A runs through all BBLs to be covered by the propagator 



8 




FIG. 8: A GE-diagram contributing to the polaron self- 
energy. It factorizes into a product of G' (t») and E(t = 
n + r c + r d + r e + Tf). 



FIG. 9: A rif-diagram contributing to the molecule self- 
energy. It factorizes into a product of r(T&) and K(r = 
t c + Td + Te). The vertical dashed lines cut (for the sake 
of better visual perception) the same r(rf,)-line. 



Gj(T',q'). We assume that q' is drawn from the distri- 
bution W-\ introduced above. In this case, the acceptance 
ratio is given by 



Prdr 



iy T (T,q)G(T',q') 
W T (r',q')G(T,q) 



n 



D v (r v ,-p v ) 



n 



J LA 



D x (t x ,Px) 



(45) 



p„ = Po + q, Pa = pa-q- (46) 

For the optimized choice of Wj, see Eq. ([55]) . 



W T (T,q)G(T',q') ^ f7(r0 
^ T (r',q')G(T,q) fi(r) 



(47) 



Diagram sign: The sign of a diagram with worms 
is somewhat arbitrary since it is not physical. The am- 
biguity is resolved by assuming that M. is always con- 
nected to I in the backward direction by an auxiliary 
unity propagator. Then, to comply with the diagram- 
matic rules, we change the configuration sign each time 
any of the following updates are accepted: (i) Reconnect, 
(ii) Open/Close updates dealing with spin-up propaga- 
tors in the forward direction, (iii) Insert/Delete in the 
molecule sector, and (iv) Redirect. For Open/Close up- 
dates dealing with spin-up propagators in the backward 
direction the sign remains the same, because here the 
sign coming from changing the number of closed spin-up 
loops is compensated by the sign in Eq. (|26p ; the same is 
also true for Insert/Delete updates in the polaron sector 
(due to our choice of the auxiliary propagator sign) . For 
precisely the same sign compensation reason, the Redi- 
rect update does change the sign despite the fact that it 
preserves the number of loops. 



E. Estimators 



Only physical diagrams with one uncovered Gj ' prop- 
agator contribute to the polaron self-energy. We will 
call them GE-diagrams. An example is shown in Fig. [5] 
[Depending on restrictions imposed on the configuration 
space (easy to implement in any scheme) physical dia- 
grams with more than one uncovered G^ propagator are 



either filtered out at the time when the contribution to 
the self-energy histogram is made, or, are not produced 
in the simulation at all.] They factorize into a product of 
G< 0) and some diagram contributing to the self-energy E. 
The utility of cyclic representation is that the uncovered 
propagator can be anywhere on the backbone. In view 
of factorization, it is easy to write the MC estimator for 
the integral (to simplify notations we omit irrelevant to 
the discussion momentum p) 



/(r)E(r)dr 



(48) 



where /(r) is some function [see, e.g., Eqs. (fill)) . (fH|]. 
Indeed, consider the estimator 



Sr 



'GT, 



fir). 



(49) 



which counts all instances of GE-diagrams with an ad- 
ditional weight /(t). Here Jgs is unity for each GE- 
diagram and zero otherwise, and r is the total duration 
in time of the E-part of the GE-diagram. The Monte 
Carlo average of this estimator is 

/>00 r-OG 

(£ I )uc oc / Gf(r')dr' /(r)E(r)dr. (50) 
Jo Jo 

Similarly, within the same scheme, we collect statistics of 
all "normalization" diagrams, see Fig. [7J using an estima- 
tor projecting to the first-order diagram with the worm, 
5 aoim . Then, 

pOO poo 
(^norm)MC « Gi / Gf '(/) dft'l Y(t)oIt. (51) 

Jo Jo 

The proportionality coefficient cancels in the ratio of the 
two averages, leading to 



/ = Gr- 



[Si 



MC 



/•CO 

/ r(r) oIt 
Jo 



(<5norm )mC 

In particular, for the optimal choice of G/v, we have 
T (Si)mc 



[ ^norm )mC 



Jo 



(52) 



(53) 



Imaginary-time integrals for the product of E(r) and 
exponentials, see Eqs. (fT3"]) and (| 14[) . is all we need to 
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determine the polaron energy and residue. For the bold- 
line generalization of the scheme described in the next 
Section we need the entire dependence of self-energy on 
time and momentum. This is achieved by differentiating 
partial contributions of GS-diagrams. For example, 



£ 



Gbhii 



(54) 



is an estimator counting contributions with r within the 
i-th imaginary-time bin of width An centered at the 
point t i . Due to linear relation between / and S(r) we 
immediately realize that (for optimized choice) 



[£s,i )mC 



'norm / MC 



An 



(55) 



In complete analogy with the polaron case, we consider 
r/\-diagrams that contain one, and only one, uncovered 
T-propagator, see Fig. [51 and use them to collect statis- 
tics for the molecule self-energy. Up to straightforward 
replacements G^ <-> T, £ <-> K all relations of this sub- 
section hold true. 



IV. SIMULATING T-MATRIX T(r,p) BY BOLD 
DIAGRAMMATIC MONTE CARLO 



FIG. 10: Diagrammatic equation for the T-matrix T(r,p). 
The arc is the vacuum spin-up propagator with the constraint 
q < kf- on its momentum q. The meaning of this equation is 
nothing but correcting the vacuum result To by subtracting 
contributions from the spin-up fermions with momenta q < 



Equation (|58|) is one of the simplest examples of prob- 
lems solvable by bold diagrammatic Monte Carlo [32| . 
We refer the reader to Ref. [32[, where the algorithm of 
solving such equations is described in detail. Here we just 
mention some specific details. We find it helpful to start 
from a good trial function for obtaining high-accuracy 
results in a short simulation time. We achieve this goal 
using the following protocol. We start the simulation by 
restricting imaginary time to be smaller than r max and 
select relatively short r max . When the result is accurate 
enough, we extrapolate it to longer times, increase r max , 
and restart the simulation with the extrapolated func- 
tion, r ext serving as the trial function, i.e. we substitute 
r = r cxt + ST to Eq. |58j) and solve for ST. If necessary, 
this procedure can be repeated several times. 



Despite relatively simple form of Eqs. ([23)) and 
tabulating the two-dimensional function r(r, p) with high 
accuracy using the inverse Laplace transform of T(uj,p) 
turns out to be a time consuming job. In this work we 
have used an alternative route based on the bold dia- 
grammatic Monte Carlo technology introduced recently 
in Ref. [§2| ■ The crucial observation is that the T-matrix 
r(r, p) can be diagrammatically related to its vacuum 
counterpart T (r,p), see Fig.[TUl with the latter is known 
analytically: 



F (t,p) 



47T 

■Tl3 / 2 



e (eF+ "- p2/4m)T ffT (T) , (56) 



were 



5t (t) = ± Ve e BT erfc(±\/^7) , (57) 



for negative/positive values of the scattering length, E — 
1/ma 2 , and erfc(a;) is the error-function. [The Fermi- 
energy and the chemical potential in Eq. (|56|) come from 
the global energy shift necessary for compliance with the 
Dyson equation shown in Fig. [TU[ ] 

With the explicit expression for the product of two vac- 
uum propagators the relation shown in Fig. [10] reads (the 
momentum argument of T(r,p) is suppressed for clarity) 

- T(r) = -r (r) + f ds f ds'T (s)T{T-s') 



'q<ki 



d Q (_ p -[(p-q) 2 +q 2 ](s'-s)/2rn 

(2tt) 3 



(58) 



V. SIMULATION RESULTS AND SERIES 
CONVERGENCE PROPERTIES 

Nearly all results in this paper were obtained by simu- 
lating diagrams built on bare one- and two-particle prop- 
agators G^ and T. We observed that the correspond- 
ing series are likely to be divergent. This, however, does 
not mean that the entire idea of calculating contributions 
from diagrams of higher and higher order and extrapolat- 
ing results to the infinite order is useless and ill-defined. 
On the contrary, it was recognized long ago that appro- 
priate re-summation techniques allow one to determine 
reliably the function standing behind the divergent series. 
Moreover, all re-summation techniques (formally, there 
are infinitely many!), if applicable, have to agree with 
each other on the final result. This important observa- 
tion vastly increases the utility of the Diag-MC technique 
we are developing here. In the next Section we demon- 
strate that making the series for £ self-consistent with 
the use of Dyson equation — bold-line technique — is an- 
other way to improve series convergence properties. 

For the resonant Fermi-polaron considered here the 
Cesaro-Riesz summation method solves the convergence 
problem. In general, for any quantity of interest — in our 
case they are polaron or molecule self-energy — one con- 
structs partial sums 



N, 



E(JV.) = D » F i 



N 



(59) 



N=l 
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FIG. 11: The molecule energy (at kF a — 1) as a function 
of the maximum diagram order TV* for different summation 
techniques: Cesaro (open squares), Riesz 5 = 2 (filled circles, 
fitted with the parabola y = -2. 6164 + 0. 28013a; + 0.01638x 2 ), 
Riesz <5 = 4 (open circles, fitted with the parabola y = 
-2.6190 + 0.61635x - 0.3515x 2 ), and Eq. flBT} (stars fitted 
with the horizontal dashed line). Reproduced from Ref. [re- 



defined as sums of all terms up to order iV* with the N- 



th order terms being multiplied by the factor F 1 



N 



In 

the limit of large 2V* and N <C iV* the multiplication fac- 
tors F approach unity while for N — > they suppress 
higher-order contributions in such a way that E(iV*) has 
a well-defined A* — > oo limit. There are infinitely many 
ways to construct multiplication factors satisfying these 
conditions. This immediately leads to an important con- 
sistency check: Final results have to be independent of 
the choice of F. In the Cesaro-Riesz summation method 
we have 



F { n* ] = [(JV.-JV + !)/#. 



(Cesaro-Riesz) . (60) 



Here 5 > is an arbitrary parameter (5=1 corresponds 
to the Cesaro method). The freedom of choosing the 
value of Riesz's exponent 5 can be used to optimize the 
convergence properties of £(A*). 

We proceed as follows. For the series truncated at or- 
der A* we first determine the polaron and molecule ener- 
gies and then study their dependence on as A, — * oo. 
In Fig. QT] we show results for the molecule energy at 
kp a = 1. Without re-summation factors the data are 
oscillating so strongly that any extrapolation to the infi- 
nite diagram order would be impossible; we consider this 
as an indication that the original series are divergent. Os- 
cillations remain pronounced for 8 = 1, but are strongly 
suppressed for larger values of 5, so that for 6 — 4 we were 
not able to resolve odd-even oscillations any more. How- 
ever, the smoothness of the curve for large 5 = 4 comes 
at the expense of increased curvature, which renders the 
extrapolation to the 1/A* — > limit more vulnerable to 
systematic errors. Empirically, we constructed a factor 
F^* ■* which leads to a faster convergence [see an example 



FIG. 12: (Color online) The polaron energy (at the uni- 
tarity point a -1 = 0) as a function of the maximum dia- 
gram order N* using Eq. (|61[). The data are fitted using 
linear -0.618 + 0.033/iV* (red) and exponential -0.6151 + 
0.026e _0,39JV * (black) functions to have an estimate of sys- 
tematic errors introduced by the extrapolation procedure. 



in Fig. (HTM 



p(N.) = C(N „) 



N, 



exp 



-N 



(A* + 1) : 



m(7V* 



(61) 



where is such that F± '> = 1. The most impor- 

tant conclusion we draw from Fig. [TT] is that in our case 
the series are subject to re-summation methods and the 
result of extrapolation is method independent. We con- 
sider small variations in the final answer due to different 
re-summation techniques and extrapolation methods as 
our systematic errors. An example is shown in Fig. [T2j 
In the next section we will present evidence that the ac- 
tual answer is closer to the upper bound of —0.615. We 
see that in the absence of additional information one has 
to allow for different ways of extrapolating the answer. 

Apart from consistency checks, one can test numeri- 
cal results against an analytic prediction for the strong 
coupling limit hp a — * corresponding to a compact 
molecule scattering off majority spins. In this limit the 
molecule energy is given by the expression 



En 



1 



e F 



2ira 
(2/3)ri 



i T (fc f o->0), (62) 



where the first term is the molecule binding energy in vac- 
uum, the second term reflects finite chemical potential of 
spin-up fermions, and the last term comes from the in- 
teraction between the composite molecule with the Fermi 
gas . The molecule-fermion s-scattering length a w 1.18a 
[35l |36| is obtained from the non-perturbative solution of 
the three-body problem. Agreement with Eq. ([62]) pro- 
vides a robust test for the entire numerical procedure of 
sampling asymptotic diagrammatic series. Our data are 
in a perfect agreement with the a « 1.18a result within 
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FIG. 13: (Color online) Polaron (black circles) and molecule 
(red triangles) energies (in units of ef) as functions of kpa. 
The dashed line on the lower panel corresponds to Eq. ((62}. 
Reproduced from Ref. 



the statistical uncertainty of the order of 5%, see the 
lower panel in Fig. 1131 

Somewhat surprising outcome is that E m is described 
by Eq. (j6"2]l very accurately all the way to the crossing 
point. This fact can be used to approximate the energy 
density functional of the superfluid polarized phase in 
the strongly imbalanced gas for kpa < 1 as that of the 
miscible dilute molecule gas coupled to spin-up fermions 
(see also Refs. HUES]) 



E = -e F n r 



1 2wa 

— ft t 

no 1 (2/3)m 1 



"1 + 



nf , (63) 



m 



where n.f, n± are densities of unpaired spin-up fermions 
and molecules, and cimm ~ 0.6a is the molecule-molecule 
scattering length [38j |. Within this approach it is found 
that the system undergoes phase separation for kpa > 

0.56 m. 

Phase separation precludes one from investigating the 
crossing point between the polaron and molecule curves 




FIG. 14: (Color online) Polaron (black circles) and molecule 
(red triangles) effective mass as functions of kpa. The vertical 
dotted line stands for (kpa) c = 1.11. The dashed line is the 
contribution from the first diagram 37]. Reproduced from 
Ref. El. 



in trapped imbalanced Fermi gases. It also rules out the 
multi-critical point on the phase diagram predicted in 
Ref. This issues, however, are not directly relevant 
to our study of properties of one spin-down particle. In 
Fig. 1131 we present polaron and molecule energies in the 
region of kpa ~ 1 where the nature of the quasiparti- 
cle state changes. The crossing point is found to be at 
(kpa)c = 1.11(2). Overall, both curves are in excellent 
agreement with the variational Monte Carlo simulations 
[lq , [l9| . There is a certain degree of accidental coinci- 
dence in the fact that the polaron self-energy is nearly ex- 
hausted by the first-order diagram considered in Ref. [37| , 
see also Fig. [T^l As we show in the next section, both 
second-order and third-order diagrams make considerable 
contributions to the answer but they happen to nearly 
compensate each other, i.e. Green's function renormal- 
ization and vertex corrections have similar amplitudes 
and opposite signs. 

The intersection of the polaron and molecule curves 
can be determined very accurately because both solutions 
are describing well-defined quasiparticles at the crossing 
point. This is because matrix elements connecting the 
two branches involve at least four particles and their 
on-resonance phase volume is zero at {kpa) c . Indeed, 
the energy, momentum, and particle number (for each 
spin direction) conservation laws dictate that polaron de- 
cays into molecule, two holes and one spin-up particle 
(molecule decays into polaron two spin-up particles and 
one hole). For this process the final-state phase volume 
gets negligibly small as compared to the energy difference 
\E p — E m \ at and in the vicinity of the crossing point. 

The data for the effective mass is presented in Fig. [TJ] 
At the crossing point the effective mass curve is discontin- 
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uous as expected for the exact crossing between two solu- 
tions. One important observation is that good agreement 
with Eq. (|62|) for molecule energy in the entire region 
kpa < 1 does not guarantee yet that the compact-boson 
approximation is accurate in the same region if proper- 
ties other than energy are addressed. In this sense, even 
good agreement with Eq. (|52^) has to be taken with a 
grain of salt since (|S2^1 assumes that the molecule mass 
is 2m and independent of kpa. The actual effective mass 
is significantly enhanced in the vicinity of (kpa) c . 



VI. NUMERIC REALIZATION OF THE 
"BOLD-LINE" TRICK 

The success of the Diag-MC method for fermi-polarons 
is, to large extent, due to small error bars we have for the 
sign-alternating sums of high-order diagrams. In gen- 
eral, it is expected that the computational complexity of 
getting small error bars for sign-alternating sums in the 
limit of large N is exponential (or factorial) in N since it 
usually scales with the configuration space volume. Any 
tricks that reduce the configuration space while keeping 
the scheme exact are worth investigating. In fact, we 
have already used one of such tricks above when we intro- 
duced r summing up all ladder diagrams. As a result, the 
diagram order was defined by the number of Gj -lines, not 
bare interaction potential vertexes, and ladder-reducible 
diagrams were excluded from the configuration space. 

In this Section we go one step further and apply an- 
other method, well-known in analytic calculations (but 
virtually never carried out analytically to high-order; for 
first-order diagrams it is known as the self-consistent 
Born approximation). It is called the "bold-line" trick. 

The relation between Gj , uf and S accounts for the in- 
finite sum of diagrams forming geometrical series. Now, 
if one-particle lines in self-energy diagrams are represent- 
ing exact Green's function (in this case they are drawn in 
bold) then many diagrams have to be excluded to avoid 
double counting. Namely, any structure which can be in- 
terpreted locally as part of the Dyson equation for G j has 
to be crossed out. Though formally the diagram order 
is still defined by the number of Gj, it is in fact repre- 
senting a whole class of diagrams (up to infinite order) in 
the original, or bare, terms. Clearly, the MC scheme is 
now self-consistently defined and potentially even finite 
number of "bold" diagrams can capture non-perturbative 
effects. 

Recently, we have demonstrated that the bold-line 
trick is compatible with Diag-MC and the corresponding 
scheme has been termed the bold diagrammatic Monte 
Carlo (BMC) [Hj]. There are two routes for implement- 
ing the bold-line technique. The first one is to arrange 
two (running in parallel) coupled Monte Carlo processes: 
one sampling the series for the self-energy in terms of ex- 
act propagators, and the other one sampling propagators 
from the Dyson equation. The latter process is essentially 
the same as the process we use for pre-calculating T, with 



an important new feature that the self-energy used in 
the sampling is permanently updated. The second route 
is specific for Dyson-type equations which allow trivial 
algebraic solution in momentum representation. In the 
present work, we use the second route. 

The implementation of the bold-line trick requires that 
we introduce an update which changes the global mo- 
mentum p of the diagram. This update applies only to 
the simplest normalization diagrams, see Fig. [7] The 
integrated weight of normalization diagrams is given by 
A(p), see Eqs. (|30[) and (|5Tj) . In the update, we select the 
new value for the global momentum from the probabil- 
ity distribution A(p) and propose new time variables for 
the spin-down and pair propagators from the optimized 
probability distributions Wr ( T i , p) an d (r 2 , p) . Since 
new variables are seeded using the exact diagram weight 
the acceptance ratio is unity. In practice, the modulus of 
the global momentum variable is defined on the discrete 
set of points. 

We start the simulation with Gj = G^\ and collect 
statistics to the momentum-time histogram of £ from 
bold-line diagrams. After a certain number of updates, 
we perform fast Fourier transform of S(r,p) to obtain 
H(u),p), calculate Gj(u>,p) using Dyson equation, which 
is then transformed back to G j (r, p) . The simulation pro- 
ceeds as before with the G±(u),p) function being recalcu- 
lated at regular time intervals to reflect additional statis- 
tics accumulated in S. Obviously, the self-consistent 
feedback present in the bold-line scheme at the beginning 
of the simulation violates the detailed balance equation 
each time the function Gj is updated. Only in the long 
simulation time limit when both £ and G do not change 
any more is the detailed balance satisfied. 

The other point which requires special care is the treat- 
ment of ladder-reducible diagrams. In the bold-line im- 
plementation we have to allow ladder diagrams back, but 
each spin-down line in the ladder-reducible structure now 
has to be understood as a difference Gj — g|°' > . Indeed, 
ladder diagrams included in T are built on bare propa- 
gators and thus have to be corrected for the difference 
between the bare and exact propagators. 

Finally, we apply the bold-line approach in the 
molecule channel as well. In fact, the scheme was de- 
signed to be identical in the one- and two-particle sec- 
tors. Now, in all self-energy diagrams we have to substi- 
tute G'j ' for Gj and T for Q, with both G and Q being 
periodically recalculated to reflect additional statistics 
statistics accumulated to the S- and X-histograms. Cor- 
respondingly, diagrams which can be locally interpreted 
as part of the Dyson equation in the molecule sector have 
to be excluded. [Ironically, this means that ladder dia- 
grams are not allowed once again with the exception for 
the first-order diagram in K which ensures that ladders 
are built on Gj, see previous paragraph.] The updates 
and acceptance ratios do not change in the bold-line rep- 
resentation, but in the optimized version the probability 
distributions Wj and Wr (and their normalization inte- 
grals) are now proportional to Gj and Q, i.e. they have 
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to be changed each time the new solutions of the Dyson 
equations are generated. 

As discussed in Ref. [23], the Monte Carlo procedure 
of solving self-consistent equations is more robust and 
has better convergence properties then standard itera- 
tions, especially for sign-alternating series. There are 
additional tools to improve the efficiency and conver- 
gence, some are self-explanatory. It definitely helps to 
start with the initial function Gj, as close as possible to 
the actual solution (the final answer should not depend 
on small variations of the initial choice). For example, 
the simulation for a given value of kpa may start with 
the final solution for the neighboring kpa point. The ini- 
tial statistics has to be discarded or "erased" according 
to some protocol. If analytic expressions in special cases 
are available, e.g. in the perturbation theory or strong 
coupling limits, they can be used to match numeric data 
and restrict the parameter range probed in the simula- 
tion. 
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FIG. 15: (Color online) Polaron energy as a function of the 
maximum diagram order iV* at the unitary point Ufcl = oo 
within the bold-line approach. Black squares show results 
when the bold-line approach was implemented only for the 
polaron propagator. When both polaron and molecule prop- 
agators in diagrams are given by G| and Q one obtains results 
shown by red circles. 

In Fig. [T5]we present data for the polaron energy at the 
unitary point calculated using the self-consistent scheme. 
As before, the diagram order is defined by the num- 
ber of spin-down propagators. To see the difference be- 
tween various approaches we first calculated E with the 
bold-line trick implemented only for spin-down propaga- 
tors (black squares) and then for both spin-down and 
molecule propagators (red circles). This plot makes it 
clear that the perturbation theory result [3 71 ] is accurate 
because corrections to spin-down propagators nearly can- 
cel vertex corrections. Figure [15] is also telling the story 
which we observe happening over and over again for other 
strongly-correlated condensed matter problems: it does 
not really make any sense to propose "better" approxi- 



mations which account only for some incomplete set of 
diagrams. For example, if in the first-order diagram V is 
replaced with Q (self-consistent Born approximation in 
the molecule sector), the answer is getting much worse! 
Moreover, using exact expressions for Gj and Q in ir- 
reducible diagrams up to third-order (!) results in an 
oscillation (the highest circle in Fig. [T5|) which forces one 
to think that the final answer is probably even further 
up. Fortunately, with the bold-line Diag-MC technique 
developed in this article we can see how different ap- 
proximations work and what their actual value is, term 
by term. 



VII. CONCLUSIONS 

The sign-problem in MC simulations is the problem of 
obtaining small error-bars for system parameters which 
allow reliable extrapolation of results to the thermody- 
namic limit. Most MC schemes are based on simulations 
of finite-size systems (of linear size L) with the configura- 
tion space volume growing exponentially/factorially with 
L 3 and inverse temperature 1/T (for quantum models). 
Since error bars grow with the configuration space vol- 
ume they are completely out of control before a mean- 
ingful extrapolation to the thermodynamic limit can be 
done. 

Computational complexity of the Diag-MC technique 
for sign-alternating series is also exponential/factorial in 
the diagram order and final results have to be extrap- 
olated to the N* — > oo limit. In this sense the tech- 
nique does not solve the sign-problem, but offers a better 
route for handling it. One important difference between 
the configuration space volume for finite-size systems and 
connected Feynman diagrams is that the latter deals with 
the thermodynamic limit directly. Moreover, the same 
diagrams describe systems of different dimensions and 
temperature. The list of advantages does not end here 
because one can employ all known analytic tools to re- 
duce the configuration space, and thus make an expo- 
nential advance towards acceptable solution of the sing- 
problem. By simulating the self-energy instead of the 
Green's function (this was not done before) the configu- 
ration space is reduced to that of G-irreducible graphs. 
Using ladder diagrams we convert the standard pertur- 
bation theory in the bare potential V into the series ex- 
pansion in terms of T. Finally, the entire scheme is made 
self-consistent by writing diagrams in terms of exact G 
and Q. Since self-consistency accounts for infinite sums 
of diagrams forming geometrical series the configuration 
space of bold-line diagrams is reduced further. All com- 
bined, the final formulation is compact enough to perform 
the TV* — > oo extrapolation reliably before error bars ex- 
plode. Strictly speaking, having convergent series is not 
a requirement because re-summation techniques are well- 
defined mathematically and their work is guaranteed by 
theorems based on properties of analytic functions. 

At the moment we do not see any obvious limitations 
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of the method described here. On the contrary, we be- 
lieve that it can be used to study a generic interacting 
many-body system, Bose or Fermi. Of course, the struc- 
ture of diagrams and the optimal strategy of applying 
analytic tools are Hamiltonian dependent and have to 
be studied case-by-case. For example, in lattice mod- 
els there is no urgent need to deal with the ultra-violet 
limit explicitly and one can proceed with the expansion 
in the bare interaction potential V; the so-called ran- 
dom phase approximation can be used to replace V with 
the screened interaction potential; the latter can be com- 
bined with ladder diagrams, etc. The bold-line trick for 
Green's functions can be implemented in any scheme. 

To summarize, we have shown that polaron type prob- 
lems can be studied numerically with high accuracy us- 
ing Diag-MC methods even when the corresponding dia- 



grammatic expansion is not sign-positive and divergent. 
Previously such series were regarded as hopeless numeri- 
cally, to such an extent that nobody was actually try- 
ing! Using Diag-MC approach we calculated energies 
and effective masses of resonant Fermi-polarons in the 
BCS-BEC crossover region and determined that the point 
where the groundstate switches from the single-particle 
(fermionic) sector to two-particle (bosonic) sector is at 
kpa = 1.11(2). This point falls inside the phase separa- 
tion region for the dilute mixture of spin-down fermions 
in the Fermi gas of spin-up particles [191 ] . 
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